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ABSTRACT 

Highly supersonic, compressible turbulence is thought to be of tantamount importance for star 
formation processes in the interstellar medium (ISM). Likewise, cosmic structure formation 
is expected to give rise to subsonic turbulence in the intergalactic medium (IGM), which may 
substantially modify the thermodynamic structure of gas in virialized dark matter halos and 
affect small-scale mixing processes in the gas. Numerical simulations have played a key role 
in characterizing the properties of astrophysical turbulence, but thus far systematic code com- 
parisons have been restricted to the supersonic regime, leaving it unclear whether subsonic 
turbulence is faithfully represented by the numerical techniques commonly employed in as- 
trophysics. Here we focus on comparing the accuracy of smoothed particle hydrodynamics 
(SPH) and our new moving-mesh technique AREPO in simulations of driven subsonic turbu- 
lence. To make contact with previous results, we also analyze simulations of transsonic and 
highly supersonic turbulence. We find that the widely employed standard formulation of SPH 
yields problematic results in the subsonic regime. Instead of building up a Kolmogorov-like 
turbulent cascade, large-scale eddies are quickly damped close to the driving scale and decay 
into small-scale velocity noise. Reduced viscosity settings improve the situation, but the shape 
of the dissipation range differs compared with expectations for a Kolmogorov cascade. In con- 
trast, our moving-mesh technique does yield power-law scaling laws for the power spectra of 
velocity, vorticity and density, consistent with expectations for fully developed isotropic tur- 
bulence. We show that large errors in SPH’s gradient estimate and the associated subsonic ve- 
locity noise are ultimately responsible for producing inaccurate results in the subsonic regime. 
In contrast, SPH’s performance is much better for supersonic turbulence, as here the flow is 
kinetically dominated and characterized by a network of strong shocks, which can be ade- 
quately captured with SPH. When compared to fixed-grid Eulerian simulations of turbulence, 
our moving-mesh approach shows qualitatively very similar results, although with somewhat 
better resolving power at the same number of cells, thanks to reduced advection errors and the 
automatic adaptivity of the AREPO code. 

Key words: hydrodynamics, shock waves, turbulence, methods: numerical 



1 INTRODUCTION 

Astrophysical gas dynamics in the interstellar and intergalactic 
medium is typically characterized by very high Reynolds numbers, 
thanks to the comparatively low gas densities encountered in these 
environments, which imply a very low physical viscosity for the 
involved gas. We may hence expect that turbulent cascades over 
large dynamic ranges are rather prevalent, provided effective driv- 
ing processes exist. Such turbulence can then be an important fea- 
ture of gas dynamics, for example providing an additional effective 
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pressure contribution, or leading to thorough small-scale mixing of 
chemical elements in the gas. 

In fact, it is believed that turbulence in the interstellar medium 
(ISM) plays a key role in the formation of ordinary stellar popu- 
lations, determining in part the initial mass function of stars, the 
lifetime of molecular clouds, and the overall efficiency of star for- 
mation (e.g. |Klessen et al.|200dl|Mac Low & Klessen|2004) . Here 
the turbulence is highly supersonic, and presumably driven primar- 
ily by supernova explosions. In addition, the strong radiative cool- 
ing processes of the ISM make its equation-of-state approximately 
isothermal, such that very strong shocks and high compression ra- 
tios are associated with the supersonic gas motions. An additional 
complexity arises from magnetic fields that are flux-frozen into the 
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gas, such that the relevant behaviour is that of isothermal, driven, 
supersonic, magnetohydrodynamic turbulence. 

Another regime where turbulence is thought to be important 
lies in cosmological structure formation, particularly in virialized 
gaseous halos, and in mildly non-linear filaments. Here small- 
scale random motions may contribute a significant fraction of the 
pressure support in group- and cluster-sized halos. For example, 
|Schuecker et al.|(2004| analysed pressure fluctuations in the Coma 
cluster, finding them to be well described by a Kolmogorov power 
spectrum with a lower limit of 10% of the total ICM pressure be- 
ing in turbulent form. Besides this direct observational evidence, 
there are also strong analytic arguments that suggest that hierarchi- 
cal mergers should be able to generate and sustain subsonic turbu- 
lence in galaxy clusters fSubramanian et al.|2006| l. 

Turbulence in the intracluster (ICM) and intergalactic medium 
(IGM) is expected to be primarily of subsonic character, allowing 
it to be approximately described as incompressible turbulence. The 
linear growth of structure driven by gravity is not expected to be 
able to generate this turbulence directly, due to the irrotational char- 
acter of the gravitational force field. However, during non-linear 



3 SPH codes) when applied to the decay of supersonic turbulence. 
They found generally somewhat higher dissipation rates in SPH 
(which can however be modified with different artificial viscosity 
parameterizations), but concluded that on scales resolved with at 
least 32 resolution elements per dimension, results were found to be 
qualitatively similar, at least as far as the decay of the Mach num- 
ber with time was concerned. The velocity power spectra, however, 
revealed that the damping on smaller scales was consistently larger 
in SPH than in the grid codes. 

In contrast, Price & Federrath (2010) claimed excellent agree- 
ment for driven supersonic Mach M = 10 turbulence between 
the SPH-code PHANTOM and the Eulerian mesh-code FLASH. In 
particular, they found a consistent Kolmogorov-like slope in the 
power spectrum of the variable p 1 ^ 3 v. Their plain velocity power 
spectra also agreed over an extended range of wave number k at 
large scales, although the SPH result eventually dipped down ear- 
lier. In the density power spectrum on the other hand, their SPH 
result yielded more power at high k, which can be interpreted as 
a welcome result of the adaptive resolution of SPH. Price & Fed- 



structure formation, curved accretion and flow shocks can intro- 


for the turbulence spectrum in SPH by Padoan et al. 


(2007 1, based 


duce vorticity into the gas. Also, the baroclinic term in the wake of 


on simulations presented by Ballesteros-Paredes et a 


. (2006), may 


mergers can act as an efficient source of vorticity, providing a large- 


have been a resolution effect only. 




scale driving of intrahalo turbulence. Already a moderate degree 


While the agreement reported by Price & Federrath|( 2010| for 



of gas turbulence may then have important consequences for the 
thermodynamic structure of quasi hydrostatic halos. For example, 
turbulence may help to create entropy cores in clusters ( jMitchell] 
|et al.|2 009 ; ZuHone 20111, and its dissipation provides for thermal 
heating. Also, the transport and the mixing of metals is affected by 
turbulent gas motions. All of these effects can modify the radia- 
tive cooling rates of intrahalo plasma, which immediately impacts 
galaxy formation at the halo centers ( [Vogelsberger et al.|201 1 [ i. In 
addition, turbulence can be expected to influence the magnetic field 
structure in virialized halo gas, which in turn affects transport pro- 
cesses in the plasma such as thermal conduction (e.g. |Parrish et al.[ 
|2012| > or physical viscosity jSijacki & Springel|2006[ l. 

A considerable difficulty for the understanding of turbulence 
lies in the comparatively limited quantitative knowledge that could 
thus far be gained from purely analytic considerations. While Kol- 
mogorov’s theory for the scaling laws of self-similar, incompress- 
ible turbulence still stands out as one of the most insightful char- 
acterizations of the physics of turbulence, an assessment of the 
accuracy of theories for turbulence, especially when constructed 
for particularly challenging cases such as compressible magneto- 
hydrodynamic turbulence, has relied to a large extent on numeri- 
cal simulations. Most simulations of astrophysical flows solve the 
Euler equations and not the Navier-Stokes equations, based on the 
realization that the residual physical viscosity is overwhelmed by 
numerical viscosity anyway at the achievable resolutions. The dy- 
namic range per dimension that can be realized in 3D turbulence 
simulations is indeed quite small, rarely exceeding a factor 10 3 at 
present. Thus only a small intertial range and comparatively low 
Reynolds numbers can be resolved directly. Nevertheless, numer- 
ous numerical studies of the properties of astrophysical turbulence 
have been carried out and have already led to significant advances 
in our understanding of this important phenomenon. 

This is especially true for the study of highly supersonic tur- 
bulence relevant for star formation, where a large body of literature 
has been accumulated, including systematic comparisons of differ- 
ent numerical techniques ( |Mac Low et al.||1998| |Kitsionas et al.| 
|2009[ |Kritsuk et al.||201 f] >. For example, |Kitsionas et al.| j2009) 
have compared different hydrodynamical codes (4 mesh codes and 



SPH and mesh codes is encouraging, it is important to keep in mind 
that this was achieved for highly supersonic turbulence. The very 
different physics of such a flow compared with subsonic turbulence 
should caution against taking it for granted that this reassuring suc- 
cess carries over to subsonic flow as well. In particular, we note that 
|Padoan et al.H2007| > observed that there appears to be a Mach num- 
ber dependence of the turbulent slope in SPH, where for M = 3 a 
slightly steeper (and hence ‘more wrong’) slope was obtained than 
for A4 = 6. This may mean that one tends to get a more accurate 
result with SPH when the Mach number is high. Indeed, in a re- 
sponse paper to our study, |PrIce| (20 1 2| pointed out that for higher 
Mach number one naturally expects a higher Reynolds number in 
SPH for given viscosity settings, which should then also allow a 
larger inertial range. 

In order to clarify these issues further, we focus in this paper 
on the topic of subsonic turbulence, which is thought to be impor- 
tant in cosmological structure formation. For example, jRyu et aT7| 
( |2008) > propose that vorticity generation along large-scale structure 
formation shocks creates significant turbulence outside filaments, 
and plays an important role in the production of intergalactic mag- 
netic fields in clusters, groups and filaments. |Dolag et al.|(|2005|l and 
|Vazza et~aT7| < |2006] > used SPH to study turbulence in galaxy clusters 
formed in cosmological simulations. They found evidence for scal- 
ing laws where the total turbulent energy in the ICM scales with 
cluster mass. Lau e t al.| (j2009) argued that the effective pressure 
associated with ICM turbulence may introduce a bias in the clus- 
ter masses inferred from hydrostatic models. [Dolag et al.[ l j2005( > 
introduced and tested a scheme for reduced viscosity in SPH, find- 
ing that this produces significantly higher levels of turbulent gas 
motions than for ordinary viscosity, reaching up to 5-30% of the 
thermal energy. In this study, they also measured a turbulent power 
spectrum for the central 500 kpc of a simulated cluster, finding a 
significantly shallower slope than expected for Kolmogorov turbu- 
lence. 

Similar results were obtained by |Valdarnini| pOTTj l, who has 
presented the most detailed study of cluster turbulence based on 
SPH thus far, including an investigation of the impact of artificial 
viscosity on the results. The study finds that a Kolmogorov-like 
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JV[ ~ 0.3 turbulence simulations 





64 3 


resolution elements 
128 3 256 3 


512 3 


SPH 


SI 


S2 


S3 


S4 


SPH (time dependent AV) 


S 1 -tav 


S2-tav 


S3-tav 




AREPO (moving-mesh) 


Al 


A2 


A3 


A4 


Fixed Cartesian grid 


FI 


F2 


F3 


F4 



Table 1. Names of primary simulation runs for A4 ~ 0.3 turbulence. 
We consider calculations with three different numerical methods. (1) SPH 
with default parameters as implemented in the GADGET code and with 
a time dependent artificial viscosity parametrization, (2) the moving-mesh 
approach of AREPO, and (3) a fixed Cartesian mesh, which is also realized 
with AREPO. 



M ~ 0.3 turbulence, variants of SPH simulations 



Run name 


resolution 


Characteristics 


S2-ngb 1 


128 3 


7V ngb = 180 smoothing neighbours 


S2-ngb2 


128 3 


Af ngb = 512 smoothing neighbours 


S3-o = 0.1 


256 3 


reduced viscosity coefficient a = 0.1 


S3-o = 0.01 


256 3 


reduced viscosity coefficient a = 0.01 


S3-ct = 0.001 


256 3 


reduced viscosity coefficient a = 0.001 


S3-ct = 0.0001 


256 3 


reduced viscosity coefficient a = 0.0001 


S3-balsara 


256 3 


enabled Balsara shear viscosity factor 


S3-tav-balsara 


256 3 


TAV + Balsara shear viscosity factor 



Table 2. Variations of the numerical parameters of our standard SPH sim- 
ulation for driven M ~ 0.3 turbulence. For resolutions of 128 3 or 256 3 , 
we carry out several simulations where either the number of SPH smooth- 
ing neighbours, or the artificial viscosity parameterization is changed, as 
indicated in the table. 

slope for the longitudinal velocity spectrum can be reached for a 
limited range of wave numbers, but the solenoidal power spec- 
trum appears to be rather strongly affected by numerical resolu- 
tion effects (in fact, the measured slopes are much steeper than 
Kolmogorov for all viscosity schemes, and this was found to be 
robust as a function of resolution over the tested range). Overall, 
|Valdamim1|201 l| l finds that turbulent motions account for a few to 
10% of the thermal energy content. 



Mesh-based studies of cluster turbulence find similar or even 
higher contributions of the kinetic energy in turbulence relative to 



the thermal energy in the ICM (Iapichino & Niemeyer||2008 , 


Vazza 


et al. 2009; Maier et al. 2009, Lau et al. 2009, Vazza et al. 


2011 


Paul et al. 201]} 


Schmidt & Federrath 2011 Iapichino et al. 


2011 


Jones et al. 2011 


. For example, Vazza et al. (201 1 ) report 5 to 30% 



based on simulations with the ENZO code. All these studies agree 
that major mergers efficiently inject turbulence, and produce a ra- 
dial trend where the importance of turbulence increases with radial 



distance. In particular, |Paul et ah] < |201 1[ ) point out that there is a 
close connection between turbulence production and virialization. 
Also, they show that the turbulence after a major merger is quite 
long-lived, still accounting for about 15% of the thermal energy af- 
ter 4 Gyrs and 5% after 10 Gyrs. |Zhu et al.| < f2010||201ll > have ex- 
amined turbulence and vorticity in the intergalactic medium, claim- 
ing that at z = 0 the IGM is in a fully turbulent state on scales less 
than about ~ 3 Mpc and that this significantly modifies structure 
formation and the gas fractions in low mass halos. 

Unlike in the case of supersonic turbulence, there is a paucity 
of systematic examinations of the accuracy of different numerical 



Mach number M ~ 1.2 M ~ 3.5 At ~ 8.4 

SPH S3-ml S3-m5 S3-ml0 

AREPO (moving-mesh) A3-ml A3-m5 A3-ml0 

Fixed Cartesian grid F3-ml F3-m5 F3-ml() 



Table 3. Simulations carried out for transsonic and supersonic driven tur- 
bulence. For the three Mach numbers examined here, M ~ 1.2, M ~~ 3.5 
and A4 ~ 8.4, we cany out simulations at a fixed nominal resolution of 
256 3 particles/cells with three different methods, SPH, a moving mesh, and 
a fixed Cartesian mesh. 





M ~ 0.3 


M ~ 1.2 — 3.5 


Ad ~ 8.4 


a 


0.014 


0.21/3.0 


12.247 


At 


0.005 


0.005 


0.005 


t S 


1 


0.5 


0.05 


^min 


6.27 


6.27 


6.27 


^max 


12.57 


12.57 


18.85 


k oc 


k~ 5/3 


jfc-s/3 


— (k - fc c ) 2 



Table 4. Summary of the parameters of the turbulent driving routine used 
in our simulations. 

techniques for representing subsonic turbulence. Only a few stud- 
ies have considered subsonic turbulence in SPH thus far l |Violeau &| 
|Issa|2007l|Monaghan|201l||Robinson & Monaghan|20 1 1 1 and we 
are not aware of a comparative analysis of three-dimensional sim- 
ulations in this regime. However, an examination of this question 
is particularly timely because serious differences between Eulerian 
mesh-codes and SPH have been reported in a number of recent pa- 
pers (e.g. Agertz et al. 2007 j t. Given also that numerical differences 
in the representation of turbulence have been suspected to signifi- 
cantly influence cooling rates of halos and therefore impact galaxy 
formation (Vogelsberger et al.|201 l||Sijacki et al.|201 l||Keres et al.| 
|2011| >, it is important to address this gap in detail. 

In this study, we therefore compare the behaviour of turbu- 
lence simulations with three different numerical methods. We use 
the smoothed particle hydrodynamics implementation of the widely 
employed GADGET code and compare it with the novel AREPO 
technique ( Springel |2010| , both using a moving mesh or a fixed 
Cartesian mesh. We primarily focus on the poorly studied subsonic 
regime, but we also perform some simulations of transsonic and 
supersonic turbulence to be able to compare our results with the re- 
cent literature, and to characterize the behaviour of our new AREPO 
code in this regime as well. In this paper, we will only examine the 
well established “standard formulation” of SPH (as described, for 
example, in ISpringel & Hemquist|2002] >, supplemented also with a 
time dependent artificial viscosity parameterization. We note how- 
ever that a number of recent works proposed extensions or modi- 
fications of classic SPH that aim to improve the accuracy of this 



method (e.g. Wadsley et al. 


2008; Price 


2008 


HeB & Springel 


2010; Read et al. 2010; Read 


& Hay field 2011 


Abel 2011). Our 



results do not necessarily extend to these new flavours of SPH, and 
it remains to be seen weather they can resolve the problems pointed 
out here. 



This paper is organized as follows. In Section [2] we outline 
our numerical techniques and describe our simulation set, as well 
as our analysis techniques. In Section[3] we present our results for 
simulations in the subsonic regime. We then turn to results for the 
transsonic and supersonic regimes in Section [4] Finally, we give a 
discussion of our findings and present our conclusions in Section[5] 
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2 METHODOLOGY 

2.1 Numerical methods 

2.1.1 Moving- and fixed-mesh simulations with AREPO 

AREPO implements a novel quasi-Lagrangian scheme for solving 
the Euler equations on an unstructured moving mesh, as described 
in detail in |SpringeIj (20 10). The mesh is defined as the Voronoi tes- 
sellation of a finite set of points that are distributed in the simulation 
volume. A finite volume approach for hydrodynamics is formulated 
on this mesh, based on a Godunov-type scheme with second-order 
accurate reconstruction and an exact Riemann solver applied to all 
interfaces for estimating hydrodynamical fluxes. 

If the mesh-generating points are kept stationary, this hydro- 
dynamical solver is equivalent to the MUSCL-Hancock second- 
order accurate scheme widely employed in many Eulerian hydrody- 
namics codes on Cartesian meshes. In fact, this equivalence can be 
made exact if the mesh-generating points are arranged on a Carte- 
sian grid, in which case the cells of AREPO also become Cartesian. 
We will carry out some of our simulations in this mode, which we 
shall refer to as “fixed-mesh”. However, the novel aspect of AREPO 
is that the mesh-generating points are allowed to move freely, with- 
out producing problematic mesh-twisting effects. In particular, the 
points can be moved with the local fluid velocity itself, thereby pro- 
ducing an adaptive, quasi-Lagrangian behaviour. In this “moving- 
mesh” mode, advection errors are greatly reduced and become in 
fact independent of the presence of a possible bulk flow, making 
the results of AREPO manifestly Galilean-invariant, unlike ordinary 
Eulerian codes. 

AREPO can additionally employ on-the-fly refinement and 
derefinement operations of its mesh, similar to adaptive mesh re- 
finement (AMR) methods. We invoke this in our moving-mesh sim- 
ulations to guarantee that the mass resolution is always approxi- 
mately constant, as in the SPH simulations that we compare with. 
To this end, cells are (de (refined if their mass deviates by more than 
a factor of two from the desired target mass resolution (which is the 
initial cell mass). We note however that such (de)refinement oper- 
ations are only rarely needed because the Lagrangian mesh motion 
already yields a nearly constant mass per cell. We also make use 
of AREPO’s mesh regularization feature, where mesh-generating 
points of highly distorted cells may receive an additional small ve- 
locity component towards the geometric center of their cell. This 
results in a more regular mesh, which reduces errors in the linear 
reconstruction step. 

We note that the AREPO code has recently been success- 
fully used in first science applications, studying first star forma- 
tion Greif et al.| 201 1) and galaxy formation ( Vogelsberge r et ah] 
|201 1 1 . There also already exist extensions to include magnetohy- 
drodynamic s l] Pa km or et al7[|201 1( >, radiative transfer ( |Petkova ~&| 
Springel 20 lT|>, as well as treatment of the full Navier-Stokes equa- 
tions l |Munoz et al.|2012| . 



2.1.2 Smoothed particle hydrodynamics 

Smoothed particle hydrodynamics (SPH) is a particle-based ap- 
proach to fluid dynamics which is popular in astronomy due to its 
geometric flexibility, automatic adaptivity, and good conservation 
properties (see e.g. Rosswog 2009, Springel 2010, for recent re- 
views). We use the simulation code GADGET-3 (last described in 
Springel, 2005} for our SPH simulations, which employs a “stan- 
dard” formulation of SPH with fully adaptive smoothing lengths 
and a simultaneous conservation of entropy and energy (Springel 



& Hernquist2002). We note however that there are many alterna- 
tive formulations of SPH, some of them quite recently proposed to 
address certain accuracy probl ems of SPH (|Pricei 2008; Rea d et ah] 
|2010| |Read & Hayfield|201 1| |HeB & Springel|2010| >. Our results 
may not necessarily apply to all of these flavours. 

In some of our simulations, we also study the influence of nu- 
merical parameters in SPH on our results, such as the number iV S ph 
of smoothing neighbours and the artificial viscosity parameteriza- 
tion. In GADGET-3, the SPH smoothing lengths hi of particles are 
adjusted such that (4-7r/3 )h^pi = N sp wfh is always fulfilled, where 
hi is the radius at which the smoothing kernel drops to zero, pi is 
the density estimate of the particle i, and m is the target mass res- 
olution (here equal to the SPH particle masses). In our default 3D 
simulations we use Ai sp h = 64 smoothing neighbours. 

The artificial viscosity is implemented as a viscous force: 

= - ^2 mjliijV iW ij , ( 1 ) 

vise j 



d Vi 

df 



where Wij is the arithmetic average of the smoothing kernels and 
Ylij parameterizes the viscous tensor. We use the following form 
l jMonaghan| 1997 ; [Springel 2005 1 ) for n, , in our default runs: 

tt a ( c » 4 ~ Cj — 3wjj) ■ Wjj 



with Wij = Vij ■ rij/\nj\ if Vij ■ rij < 0, otherwise Wij = 0. For 
this definition of Wij, the artificial viscosity is always repulsive, and 
is non-zero only if a pair of particles approaches each other, imply- 
ing that the entropy produced by the viscosity is positive definite. 

One general problem of artificial viscosity parameterizations 
is that they may introduce spurious viscosity also outside of shocks, 
in regions where it should in principle not be needed (e.g. Cullen 
& Dehnen 2010). This can be a significant problem in shear flows, 
where this effect may lead to unwanted angular momentum trans- 
port. To suppress the artificial viscosity in regions of strong shear, 
|Balsara| ( | 1995| > proposed a simple viscosity limiter in the form of an 
additional multiplicative factor ( fi + fj)/2 for the viscous tensor, 
defined as 



fi = 



\\7-v\i 

|V ■ v\i + |V x v\i ' 



(3) 



This limiter is often used in cosmological SPH simulations and also 
available in the GADGET code. In our default simulations with fixed 
a, we have refrained from enabling it, but we have also run com- 
parison simulations where it is used, as discussed in our results 
section. 

In addition, we consider a so-called time variable artificial vis- 
cosity, as first proposed by |Morris| ( |1997j l. The idea is here to try to 
reduce the viscosity in regions away from shocks such that it is ap- 
plied in a more targeted fashion only where it is really needed. We 
employ the implementation of |Dolag et al.| < [2005} ) in the GADGET 
code, where a is replaced with an individual parameter «i(f) for 
each particle: 



dtti 

df 



Ot-i C^min 
T 



+ Si. 



(4) 



The time evolution is controlled by a source term Si that ramps 
up the viscosity quickly when a fast compression is detected, and 
a decay function that makes the viscosity decline again in smooth 
regions of the flow to a small minimum value a m i n over the decay 
time r. We note that a near-identical formulation has also been used 
by [Pnce|p0l2) . 
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2.2 Turbulent driving 



In this work, we consider isothermal gas in which turbulence is in- 
duced through an external stochastic forcing on large scales. The 
condition of isothermality is not crucial for our study of subsonic 
turbulence, but it conveniently prevents that the turbulent kinetic 
energy dissipated in the flow leads to a gradual increase of the 
pressure in the gas with time. Instead, the dissipated energy is sim- 
ply lost from the isothermal system, so that a statistically quasi- 
stationary state of developed turbulence can be reached after some 
time, where on average the energy injected on large scales is lost 
on smaller scales by dissipation. 

Our method for calculating the acceleration field follows 
Schmidt et al. |2006jl; FederraAetalJ 
(2010| l and |Price & Federrath| ( |201Cty . 



closely the procedure used in 



(2008 2009l>; |Federrath et al. 



In particular, the acceleration field is setup in Fourier space and 
only contains power in a small range of low frequency modes be- 
tween fc m in = 6.27 and fc max = 12.57. The relative amplitude of 
the forcing modes over this small range is varied as P(fc) oc k~ 5 ' 3 . 
Except in our run at Af ~ 8.4, P(k) is a paraboloid centred around 
(fc m in + fcmax)/2 with fc m in = 6.27 and fc max = 18.85. The phases 
of the Fourier modes are drawn from an Ornstein-Uhlenbeck pro- 
cess and are periodically updated after a time interval At. The cor- 
responding random sequence is given by 



x t 



= f Xt-At + £?V (1 - f 2 ) ; 



(5) 



where / is a decay factor given by / = exp(— Af/f s ), with t s 
being the correlation length. z n is a Gaussian random variable and 
a is the variance of the Ornstein-Uhlenbeck process. The resulting 
sequences have zero mean, (x t ) = 0 , and their correlations are 
given by {xt Xt+At) = er 2 /. The frequent but correlated changes 
of the acceleration field as a function of time result in a smoothly 
varying turbulent driving field. 

We use a purely solenoidal driving in this study, which can be 
obtained by projecting out the compressive part of the acceleration 
field through a Helmholtz decomposition in fc-space. The projec- 
tion operator is given by 

a{k)i = (da - a 0 (k)j ( 6 ) 

in Fourier space. We note that solenoidal driving appears partic- 
ularly appropriate for subsonic turbulence. Compressive modes 
would only cause additional sound waves and only start to cou- 
ple to smaller modes once non-linear steepening of these acoustic 
waves becomes important. In any case, if a compressive component 
was added, we would expect a somewhat broader density PDF for 
a given Mach number (Federrath et al. 2008}. 

Finally, the acceleration field due to the driving mechanism is 
calculated in position space at each particle or cell position directly 
as a sum over the small number of non-zero Fourier modes, a pro- 
cedure which is free of any resolution limitations. This field is then 
introduced as an additional source term in the Euler equations, in 
the same way as an external gravitational field is normally coupled 
to gas dynamics. In our time-integration scheme, the external ac- 
celeration is added in two half-steps at the beginning and end of 
each timestep, producing a leap-frog type integration scheme. 



2.3 Simulation set 

We consider periodic boxes of unit length on a side, filled with 
isothermal gas at unit mean density and unit sound speed (p = 
c 3 = 1). The Euler equations are scale-independent, so that once 



quasi-stationary turbulence has developed the only characteristic 
parameter remaining is the mean Mach number Ai , which we de- 
fine as the mass-weighted rms-velocity relative to the sound speed. 
The amplitude of the driving field determines Al, and can in prin- 
ciple be freely adjusted to reach the desired strength of turbulence. 
We note that we have here not attempted to subtract the mean gas 
velocity in our simulations even though a respectable amount of 
bulk motion can be generated in the supersonic regime as a result 
of our driving. In fact, the kinetic energy of the bulk motion can be- 
come at times nearly as large as the kinetic energy of the irreducible 
smaller scale motions, and hence matters when measuring the mean 
Mach number in terms of the total kinetic energy. The bulk motion 
is negligible in the subsonic regime. As the bulk motion however 
only affects the DC mode in Fourier space, it has no influence on 
our power spectrum measurements. The Mach numbers we report 
are corrected for the bulk motion and do not include it. 

For technical reasons having to do with our measurement tech- 
nique for dissipation discussed below, we actually do not use an ex- 
act isothermal equation of state, but rather one with an adiabatic in- 
dex of 7 = 1 . 001 , combined with enforcing the entropy of the gas 
to stay at the initial value after completion of each timestep. Specif- 
ically, for a prescribed initial mean density p, and an (isothermal) 
sound speed c s , we initialize the gas with a specific entropy 

A=clp'-\ (7) 

hence the pressure is given by 

P, = PCs ( ? V (8) 

for a cell/particle of density p, . The specific internal energy per unit 
mass, Ui, can be calculated from the specific entropy as 

'y — 1 

Ui = A/t—. (9) 

(7-1) 

In our ‘quasi-isothermal’ simulations, the entropy of the gas is reset 
to A after each timestep, so that the pressure of a particle/cell is 
always given by equation §. All our simulations are started from 
a regular Cartesian grid of particles with initially zero velocities, 
and use a global time step for all particles/cells. 

Each of our primary subsonic simulations was performed with 
three different numerical methods: SPH as implemented in the 
GADGET-3 code, moving-mesh hydrodynamics using AREPO, or 
fixed-mesh Eulerian hydrodynamics also based on AREPO but with 
a stationary Cartesian mesh. The resolution of our simulations 
ranges form 64 3 to 512 3 particles or cells, respectively. Table [T] 
gives an overview of these simulations, and lists their most impor- 
tant parameters. Our naming convention is such that a leading cap- 
ital letter indicates the type of a simulation (‘S' for SPH, ‘A’ for 
moving-mesh with AREPO, and ‘F’ for a fixed mesh), followed by 
a digit that indicates the resolution level (‘F for 64 3 , ‘2’ for 128 3 , 
etc.). When we compare different simulation techniques or differ- 
ent numerical resolutions, we always do this at identical driving 
amplitude, such that any difference that is seen arises from the hy- 
drodynamics alone. In particular, all of our subsonic simulations 
(where Ai ~ 0.3) use an identical turbulent forcing field. 

In our default SPH simulations, we consider two different ar- 
tificial viscosity formulations, one with a viscosity strength param- 
eter equal to a = 1 , the other with a time- variable viscosity that 
is individual for each particle (both with and without a shear vis- 
cosity limiter). These choices approximately bracket the range of 
viscosity settings that are in use in production calculations in cos- 
mology. In order to examine the dependence of our results on nu- 
merical parameters of SPH in more detail, we have additionally 
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performed a set of SPH simulations with further variations in the ar- 
tificial viscosity parameters (a = 1.0 with a shear viscosity limiter, 
a = 0.1, a = 0.01, a = 0.001, and a = 0.0001). Furthermore, 
we have also studied the influence of the number of SPH smooth- 
ing neighbours. Our default value for 7V ng b is 64, but we also used 
IVngb = 180 and iVngb = 512. The corresponding simulation runs 
and their symbolic names are summarized in Table[2] 



2.4 Measuring dissipation 

The classical theory of Kolmogorov for incompressible turbulence 
conjectures that energy is injected on large scales and then cascades 
down to eddies of ever small size, until dissipation on very small 
scales eventually occurs. In this picture, large scale gas motions 
in the resulting turbulent cascade do not dissipate energy in any 
significant way, instead the energy is essentially transported in a 
conservative fashion over the inertial range down to the dissipation 
scale. In the analysis of numerical turbulence simulations it is stan- 
dard procedure to examine how the kinetic energy is distributed as 
a function of scale, which is usually done in terms of the velocity 
power spectrum. We suggest here that it is also interesting to try 
to directly measure the energy dissipation as a function of scale, 
as this provides interesting complementary information about the 
dissipative properties of a numerical scheme. 

To this end, we define dissipation as the irreversible conver- 
sion of kinetic energy into heat. Because we use 7 = 1.001 in 
our ideal gas equation-of-state, dissipation manifests itself as an in- 
crease in the specific entropy of a cell or a particle. To measure this 
quantity, we compare Ai after completion of every timestep with 
the value A, afterwards resetting A, to A. The extracted thermal 
energy, A Ei = rrii(Ai — A)pJ~ /(7 — 1), is then the dissipated 
energy, which we assume to leave the system in concordance with 
the quasi-isothermal conditions we impose. 

In SPH, Ai only changes due to the artificial viscosity and 
is easily obtained from the work done against the artificial viscos- 
ity forces. The A Ei measured for a SPH particle is always positive 
definite in this case. In our mesh-code AREPO, exact energy-, mass- 
and momentum-conservation is given in every hydrodynamic step. 
To measure the dissipative increase of entropy of a cell, we ad- 
ditionally advect the entropy in the system in a conservation fash- 
ion, as described by Springel ( |2010j >. The dissipative energy change 
A Ei of a cell can then be estimated in the above fashion, with the 
caveat that the energy A Ei is not guaranteed to be positive definite 
for all cells due to discretization errors. However, in this case a lo- 
cal average over a group of cells will still give a faithful estimate 
of the total energy that is lost, due to the manifestly conservative 
properties of the mesh-based evolution of the fluid. We shall use 
the A Ei values for a Fourier analysis of the spatial scales on which 
most dissipation occurs, and for cross-checking whether the total 
dissipated energy balances the total injected energy when steady- 
state turbulence is reached. 



2.5 Power spectrum measurements 

The 3D power spectrum of a scalar or vector field w is defined as 
the Fourier transform of the two-point correlation function 

C w (l)= (w(x + l)w{x)) m . ( 10 ) 




Figure 1 . Results for different approaches to measure the velocity power 
spectrum for one of our t i 4 ' SPH simulations. The green lines with different 
line styles show our nearest point sampling, with sampling resolutions 32 s , 
64 3 . 128 3 , 256 4 and 512 3 . The results for 128 3 (our default grid size at this 
resolution) and higher are virtually identical up to the Nyquist frequency 
of the run. The red line shows a measurement when the velocity field is 
calculated by SPH-smoothing for a 128 3 mesh. This method (advocated by 
Pric el20l2) suppresses small-scale velocity noise but also removes kinetic 
energy associated with particle motions on these scales. 

Thus 

£ 4 fc) = 7 -L- [ C w (l)eM~ikl)d 3 l ( 11 ) 

(Zn)° J v 

= (t) i^i 2 ’ (12) 

where w is the Fourier transform of w. In order to numerically es- 
timate w, the field w is usually represented in discretized form on 
a Cartesian grid, allowing an efficient measurement of the power 
spectrum through discrete Fourier transforms. 

We here use a nearest neighbour sampling of the intrinsic hy- 
drodynamical quantities of the particle/cell data of our simulations 
to do this. The value of the desired quantity w at each grid cell 
is hence obtained as the value of the particle or cell closest to the 
centre of a grid cell. This simple approach aims to maximize the in- 
formation content extracted from the simulations, but risks to suf- 
fer from power aliasing effects if the employed Fourier grid is too 
small. We have however checked that our default Fourier mesh size 
we employ for our power spectrum measurements (based on a grid 
twice as fine as the resolution of the simulation at hand) is fine 
enough to make such effects negligible. 

Another important aspect of our power spectrum measurement 
is that it faithfully represents the total kinetic energy of the parti- 
cle/cell set. This should be given as the integral over the measured 
power spectrum. According to Parseval’s theorem, this integral is 
equal to the variance of the velocity field that is used to measure 
the spectrum through a discrete Fourier transform. As our veloc- 
ity field definition creates a fair sample of the particle/cell velocity 
values, we obtain an accurate accounting of the total kinetic en- 
ergy in the flow. In contrast, smoothing the SPH velocity field via 
kernel interpolation as advocated by |Price| ( |20 1 2] > removes kinetic 
energy on small scales and can cause a significant error in the to- 
tal energy represented by the power spectrum measurement. This is 
explicitly demonstrated in Figure |T| where we show a power spec- 
trum measurement for one of our 64 3 SPH simulations carried out 
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with different techniques. The green lines with different line styles 
show our nearest point sampling, with resolutions 32 3 , 64 3 , 128 3 , 
256 4 and 512 3 . The results for 128 3 and higher are virtually iden- 
tical up to the Nyquist frequency of the run and accurately reflect 
the total kinetic energy of the particle set. On the other hand, if we 
SPH-smooth the velocity field on a 128 3 mesh, we obtain the red 
solid line. While this method suppresses small-scale velocity noise, 
it also reduces the kinetic energy in the field by ~ 3.6% in this 
case. However, if for example our Sl-tav simulation is considered 
the suppressed energy amounts to about ~ 13% of the total kinetic 
energy due to the higher power on smaller scales in that simula- 
tion. We note that we calculate our power spectrum measurements 
on-the-fly while the turbulence simulations are run, allowing us to 
reach a very fine temporal resolution for the evolution of the power 
spectrum of the different quantities we examine. 

Assuming an isotropic statistical distribution, the ID power 
spectrum of the quantity w can then be obtained by angular aver- 
aging E w ( k ) over shells in fc-space, as 

E w (k) = 4nk 2 {E w (k)) (13) 

where k = \k\. We employ fine logarithmic bins in k for deter- 
mining the mean power per mode at a certain k, with bins com- 
bined as needed such that a minimum number of modes per fc-bin 
is obtained. The normalization of the ID power spectrum is cho- 
sen such that the integral over the power spectrum gives always the 
total power, i.e. 

1 Ar_1 

o 2 = / E w (k)dk = — ^ \wijk\ 2 . (14) 

J j,k,l = 0 

We note that we usually plot the quantity kE w (k ) in our 
power spectrum plots as a function of log k, instead of using E w ( k ) 
directly. This has the advantage that a horizontal line in such a plot 
corresponds to constant total power per logarithmic decade, one 
with a positive slope means that small scales dominate, whereas a 
negative slope indicates that the total power in the examined quan- 
tity is dominated by large scales. 

Finally, we want to clarify how we measure the power spec- 
trum of the energy dissipation rate, which requires a special treat- 
ment. In order to allow for a direct comparison with the kinetic 
energy power spectrum obtained from w = v, we in principle want 
to set w = v A E, where A E is the energy dissipation rate mea- 
sured as described above. However, since the measured dissipation 
rate can sometimes exhibit negative values in the case of the mesh 
code, this procedure needs to be modified. We thus compute the 
power spectrum for w+ = y/AE+ and = \/AEZ separately, 
where A E + = min(A E, 0) and A E- = max(A£, 0) are the 
positive and negative parts of the measured dissipation field. Fi- 
nally, the dissipation power spectrum is then estimated as 

E diss [k ) = E w+ (fc) - E w _ ( k ). (15) 

Note that the fc-integral of this quantity is equal to the total energy 
dissipation rate. 



3 SUBSONIC TURBULENCE 
3.1 Global characteristics 

In Figure [2] we show the time evolution of the rms Mach number 
for our runs of subsonic turbulence at a resolution of 256 3 (S3, S3- 
tav, A3, and F3). After an initial ramp up of the turbulent energy, 
a quasi-stationary state is established, starting at time t ~ 5 — 10. 




Figure 2. Mean Mach-number evolution as a function of time for our sub- 
sonic turbulence simulations. Here M is defined as the mass-weighted rms 
velocity in units of the sound speed. The different lines give results for SPH 
(green and yellow), AREPO (blue), and a fixed-mesh (magenta), at a res- 
olution of 256 3 particles/cells (runs at different resolutions give extremely 
similar results). We see that a quasi-stationary state is reached after time 
t ~ 10, but the total kinetic energy in the a = 1.0 SPH case (green) 
is somewhat smaller than in the two mesh codes and in the SPH run with 
time- variable artificial (tav) viscosity (yellow). 




Figure 3. Time evolution of the total cumulative injected energy (dashed 
lines) and the total dissipated energy (solid lines), for different simulation 
runs, as labeled. At any given point in time, the difference between the 
injected energy and the dissipated energy is the kinetic energy stored in the 
simulation box. The time average energy dissipation rate per unit mass is 
e ~ 0.016. 



There are however still substantial intermittent fluctuations in the 
global rms Mach number, making it clear that averaging over ex- 
tended periods of time is required to obtain truly stable results for 
the statistical properties of the turbulent fluid state, especially on 
large scales. We note that runs carried out with different numeri- 
cal resolutions give extremely similar results to the ones shown in 
Fig. 0 Interestingly, the time evolutions of the moving-mesh and 
the fixed-mesh results agree very well with each other, but the ter- 
minal Mach number reached by SPH is significantly lower, espe- 
cially in the run with a = 1.0. This is despite the fact that the 
driving field imposes exactly the same accelerations in all the sim- 



© 0000 RAS, MNRAS 000, 000-000 



8 A. Bauer and V. Springel 





SPH 



SPH 



SPH 




Figure 4. Visual comparison of the turbulent velocity field (top row), the density field (middle row) and the enstrophy V x v\ 2 (bottom row) in quasi-stationary 
turbulence with ,V1 ~ 0.3, simulated with different numerical techniques. Shown are thin slices through the middle of the periodic simulation box at the final 
time t = 25.6. From left to right, we show our moving grid result (A3), an equivalent calculation on a static mesh (F3), and two SPH calculations, one with a 
fixed a = 1.0 viscosity (S3) and the other with time-variable viscosity (S3-tav), as labeled. 



illations. The smaller overall kinetic energy achieved in the SPH 
run with a = 1.0 is a result of viscous damping of large-scales 
modes at or close to the driving scale. This effect is greatly reduced 
but not completely eliminated with the time-variable viscosity pa- 
rameterization. 

We show the cumulative injected and dissipated energy as a 
function of time in Figure [3] for the same simulations. Note that 
the difference between these two quantities is exactly the kinetic 
energy stored in the gas at the corresponding time. Interestingly, 
the mesh-based simulations do hardly dissipate any energy until 
t = 5, in contrast to the SPH simulations which show signs of 
energy dissipation right from the start. This is consistent with the 
impression from Figure |2]that it is harder in SPH than in the mesh- 
code to set the largest eddies into motion. At around t ~ 13, the 
total cumulative dissipated energies begin to be rather similar for all 
three methods, but the total injected energy of the SPH simulations 
still lags behind the mesh-based runs. This is simply because the 
lower velocities in SPH reduce the average value of (itadriv), where 
ddriv is the acceleration due to the external driving field. As a result. 



the kinetic energy in the SPH runs never fully manage to close the 
gap to the mesh-based calculations. 

3.2 Visualizations of the turbulent velocity, kinetic energy 
density and vorticity fields 

To better understand the systematic differences between the simula- 
tion techniques, it is instructive to consider maps of fluid quantities 
in slices through the simulation cube. To construct them, we use 
nearest neighbour sampling at the coordinates of a two-dimensional 
grid with 512 2 pixels, twice finer than the nominal resolution of 
the simulations examined here. Each pixel will hence show the 
value of the closest mesh cell or SPH particle, respectively. This 
directly reflects the individual fluid elements used in the discretiza- 
tion schemes of the two numerical techniques, highlighting mesh 
or sampling artefacts if they exist, as well as discretization noise. 

In the top row of Figure [4] we show slices of the velocity 
field at the final time t = 25.6 of our subsonic turbulence simula- 
tions, comparing the moving-mesh calculation with AREPO to the 
one with a fixed Cartesian mesh, and to our SPH simulations with 
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Figure 5. Velocity power spectrum of SPH and AREPO, compared at a 
resolution of 256 3 . The thin grey line shows the slope expected for Kol- 
mogorov’s theory of incompressible turbulence. The power in SPH falls 
much more rapidly than expected for fully developed turbulence. On small 
scales, the power rises again quite strongly up to the Nyquist frequency. 
This is small-scale velocity noise characteristic of SPH. The vertical green 
dotted line indicates the scale 27r//i max , where /i max is the maximum SPH 
smoothing length among all the particles. 



a = 1.0 and time-variable viscosity. We can see that the slices 
of the moving mesh (top left panel) and the fixed mesh calcula- 
tions (top middle panel) appear qualitatively very similar, featuring 
both large-scale coherent motions and many irregular small-scale 
features. The moving-mesh result appears slightly less smooth and 
shows more small-scale features, but based on these images alone 
it would be hard to decide whether this is due to a higher effective 
resolution or due to the more irregular shaped cells in the Voronoi 
case, which may induce some aliasing effects in the pixelized map. 

In contrast, the nearest neighbour SPH results for the velocity 
field, shown in the two panels on the top right, look dramatically 
different. Here the corresponding velocity fields do not contain the 
small-scale velocity features present in the mesh-based calcula- 
tions, particularly in the a = 1.0 run, suggesting that an equally 
well developed turbulent cascade has not really formed in the SPH 
simulation. 

This impression is compounded by slices through the density 
and enstrophy fields (| V x w| 2 ), shown in the middle and bottom 
row of Figure [4] respectively. While the mesh-based calculations 
exhibit a delicate mix of fine structures in the kinetic energy density 
both on large and small scales, the SPH a = 1.0 simulation shows 
only some large-scale motions, presumably reflecting primarily the 
driving field. While the S3-tav run shows more structures, there is 
still a paucity of smaller flow features. Similarly, whereas the en- 
strophy slices reveal a granular structure in the vorticity field that is 
dominated by small structures, these are essentially completely ab- 
sent in the SPH calculations. The utilization of a time dependent ar- 
tificial viscosity scheme improves the SPH result noticeably. A fur- 
ther improvement might be achieved through more advanced vis- 
cosity switches, like those suggested by|Cullen & Dehnen| j2010|) 
or |Read & HayfieldlpoTT) . 

3.3 Velocity power spectra of subsonic turbulence 

A more quantitative analysis of this difference is obtained by con- 
sidering velocity power spectra of these four different simulation 



techniques. In Figure [5] we compare power spectra of the kinetic 
energy for our runs at 256 3 resolution, averaged over an extended 
period of time from t = 10 to t = 25.6, after a quasi-stationary 
turbulent state has been established. 

The results confirm the impression inferred from the previ- 
ous subsection. There is a rather striking difference between the 
mesh-based simulations and our SPH calculations. The kinetic en- 
ergy in SPH is already drained at rather large scales, such that an 
extended energy cascade is not formed. The self-similar turbulent 
power spectrum expected based on Kolmogorov’s theory for in- 
compressible turbulence ( E v (k ) oc fc _5//3 ) is indicated as a thin 
grey power law in the figure - this line has a different slope com- 
pared with the rapid decline of the power spectrum observed in 
SPH. 

In contrast, the mesh-based simulations show a slope simi- 
lar to the expected Kolmogorov law, at least on very large scales. 
Fits to the power-law region of the velocity power spectrum give 
slopes of —1.64 and —1.68 for the moving mesh and the fixed 
mesh, respectively, whereas SPH shows a wrong slope of —4.14 
in the case of a constant artificial viscosity and —2.1 in the case 
of a time-dependent artificial viscosity. There is even an excess of 
power in the mesh-based results above the expected continuation of 
the Kolmogorov slope, before the velocity power spectrum eventu- 
ally starts to rapidly decline on scales somewhat larger than the 
Nyquist frequency corresponding to the nominal spatial resolution. 
This bump in power is a manifestation of the so-called bottleneck 
effect, which is commonly encountered close to the resolution limit 
in mesh-based studies of turbulence and is also seen in experiments 
(e.g. |Meyers & Meneveau|2008| and references therein). The nu- 
merical bottleneck effect is similar to the physical bottleneck ef- 
fect and considerably complicates attempts to robustly measure the 
true slope of the inertial range of turbulence, as this requires the 
use of extremely high resolution (2048 3 and beyond), such that the 
bottleneck bump moves to sufficiently small scales. We note that 
the moving-mesh and fixed-mesh calculations agree well with each 
other up to quite high k, where eventually some small differences 
arise, caused by the different mesh geometries and truncation errors 
in the two schemes. 

Another interesting feature of the SPH velocity power spec- 
trum is that there is actually a minimum at some intermediate 
scale, followed by a strong rises towards still smaller scales. The 
minimum occurs on scales that should formally still be well re- 
solved, because these scales are considerably larger than the max- 
imum SPH smoothing length fi max among all the particles (indi- 
cated as the vertical green dotted line at 27r//?. max in Fig. [5}. Nev- 
ertheless, already on this comparatively large scale, the power starts 
to increase again. This is a result of the high small-scale subsonic 
velocity noise present in SPH that we already witnessed in Fig- 
ure [4] Even when a kernel-smoothed velocity field is considered 
instead, and the power spectrum is calculated for this smoothed 
quantity, there is a considerable small-scale bump left, as shown 
by the dashed green line which shows the power spectrum of the 
SPH-smoothed velocity field. On large scales, the behaviour of this 
field is the same as for the nearest neighbour interpolated one, as 
expected. 

In Figure[6] we show a resolution study for the subsonic veloc- 
ity power spectra of our AREPO and SPH runs, ranging from 64 3 to 
512 3 particles/cells. The SPH simulations seem to converge to each 
other only on the largest scales. Even with a resolution as high as 
512 3 particles, there is build up of an extended inertial range with 
the expected energy cascade. We only see that with improving res- 
olution there is a slight shift towards smaller scales of the rapid 
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Figure 6. Convergence study for the velocity power spectrum of M ~ 0.3 
subsonic turbulence. The panel on top shows results for AREPO, from a 
resolution of 64 3 to 512 3 cells. The panel on the bottom gives the corre- 
sponding results for SPH. However, even at a high resolution as high 512 3 
particles, no extended inertial range of turbulence can be identified in SPH. 
The thin grey lines show the power-law expected for Kolmogorov’s theory. 



decline of the power spectrum. Also, the minimum of the power 
spectrum is reached at progressively smaller scales, but the over- 
all shape of the velocity power spectrum does not improve signifi- 
cantly, and the small-scale noise bump remains present. 

For the simulations with AREPO, we observe that the bottle- 
neck effect moves to smaller scales with improving resolution. This 
is expected, as this effect should be tied to the numerical dissipation 
occurring on scales close to the resolution limit. As the bottleneck 
moves towards smaller scales, a larger inertial range with a self- 
similar power-law region is established on large scales. We note 
that the rise of the power in the moving mesh-code on very small 
scales, at around the Nyquist frequency, is due to noise and alias- 
ing effects at the spatial resolution limit that is reached here, which 
is qualitatively a very different effect from the small-scale velocity 
noise that sets in in SPH on much larger scales. 

The computational cost of one of our moving-mesh turbulence 
simulations with AREPO in the sub-sonic regime is nearly a factor 
of 4 higher than a corresponding SPH simulation with GADGET-3 
at the same number of resolution elements. Compared to a cor- 
responding fixed-mesh calculation, the moving-mesh simulation 
is about 5 times slower. This difference arises mainly due to the 
costly Voronoi mesh construction, and in the case of moving ver- 




Figure 7. Dissipation power spectra for AREPO and SPH runs at differ- 
ent resolutions, compared to the corresponding shape of the velocity power 
spectrum at 256 3 resolution (dashed lines). For the mesh-code, the dissi- 
pation actually peaks on scales where the power spectrum starts to deviate 
from Kolmogorov’s self-similar scaling. In contrast, SPH shows very strong 
dissipation already on larger scales, preventing the build-up of a turbulent 
cascade. In addition, the dissipation is also strong on small scales, close 
to the resolution limit, where the small-scale noise developing in SPH is 
constantly damped away. 



sus fixed mesh, additionally due to the roughly twice as many Rie- 
mann problems that need to be solved for the unstructured Voronoi 
mesh compared with a Cartesian grid. However, these differences 
in run time of order unity are implementation dependent and ul- 
timately of limited importance. What should really be considered 
is the computational cost to reach a desired level of accuracy. For 
example, as our moving mesh run with 128 3 cells (A2) produces 
a model for the Kolmogorov spectrum at least as good as the SPH 
run with 256 3 particles (S3), one may state that AREPO is actu- 
ally about 4 times as efficient as SPH when comparing these two 
runs. But we caution that such statements about the relative effi- 
ciency of different schemes are in general resolution and problem 
dependent. For example, as we shall argue below (subsection |3.6| l, 
we expect that the effiency gain of AREPO relative to SPH at fixed 
accuracy actually grows with Reynolds number. Also, the relative 
cost of moving-mesh vs. fixed-mesh depends on the Mach number 
of the turbulence, because for highly supersonic flows consider- 
ably smaller timesteps are needed for a fixed mesh compared with a 
moving mesh. Finally, we remark that in astrophysical applications 
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that include self-gravity the run time difference between AREPO 
and GADGET at the same number of resolution elements shrinks 
considerably, simply because of the comparatively high cost of the 
tree-based gravity calculation. 

3.4 Dissipation power spectra 

In Figure[7] we show power spectra for the energy dissipation rate, 
measured as described in Sections |2.4| and |2.5| In the top panel, 
we show results for the simulations A1 to A3 with resolutions 64 3 
to 256 3 , averaged over the same period of time as in our veloc- 
ity power spectrum plots. For comparison, we also plot the kinetic 
energy power spectrum as dashed lines, in order to allow a compar- 
ison of the shapes of the different curves. Interestingly, the AREPO 
simulations show a peak in dissipation right at the scales where 
the velocity power spectrum begins to rapidly fall. While there is 
also some residual dissipation at very large scales (which becomes 
smaller with better resolution), this is more than an order of mag- 
nitude lower than the energy drained around the scales where the 
dissipation measurement peaks. The result is hence consistent with 
an interpretation where only negligible dissipation occurs on large 
scales, with all the energy dissipated on some smaller dissipation 
scale, which in our case is related to the numerical resolution limit. 
Such a scenario is consistent with the theoretical assumptions that 
enter Kolmogorov’s theory of self-similar scaling. 

In the bottom panel of Figure [7] we show the corresponding 
SPH results. Here a very different shape of the dissipation power 
spectrum is found. There is a peak already on very large scales, 
close to the driving scale. The amplitude of the dissipation lies 
considerably higher on these scales range than in the mesh-code, 
and shows only a weak dependence on numerical resolution. This 
explains why there is not much energy left to be fed into a turbu- 
lent cascade that could transport it conservatively towards smaller 
scales. Interestingly, there is however a second extended maximum 
of the SPH dissipation power spectrum on very small scales, co- 
inciding with the location of the small-scale bump in the velocity 
power spectrum. This is apparently related to viscous dissipation of 
some of the small-scale velocity noise in SPH. 

3.5 Dependence on SPH parameter settings 

Given the sobering results we have thus far obtained for subsonic 
turbulence in SPH, it is an important question whether this outcome 
can be significantly improved with different parameter choices for 
the method. The primary numerical parameters that may strongly 
affect the SPH results are the number of smoothing neighbours, and 
the artificial viscosity parameterization. In fact, these are the only 
aspects that can be changed easily without reverting to an entirely 
different formulation of SPH, or a substantially different method 
for particle hydrodynamics. 

An increase in the number of smoothing neighbours should 
reduce the noise in SPH kernel estimates. In fact, it has been ar- 
gued that convergence of SPH requires a simultaneous increase 
both of the number of simulation particles and a (slower) increase 
of the number of smoothing neighbours l |Rasio |2000| |Robinson &| 
Monaghan 2011). Unfortunately, in practice the clumping instabil- 
ity present for the normal SPH kernel shape counteracts attempts to 
improve the SPH estimates through a drastic increase of the number 
of smoothing neighbours (but see Read et al.|2010} . Regardless, we 
have examined whether an increase of the number of neighbours to 
Angb = 180 or even Ai ng b = 512 improves our results. To this end 
we have repeated our S2 simulation with these settings. 





Figure 8. Effects of a larger number of SPH smoothing neighbours. The 
panel on top gives results for the velocity power spectrum when the number 
of neighbours is increased from our default of 64 to 180, and finally to 512. 
Formally, the later run with 128^ particles has the same mass and spatial 
resolution as our SI run with 64 ’’ particles. The lines shown in grey cor- 
respond to an early time, when the turbulence is not yet fully established. 
Here greater differences between the results can be seen, with the larger 
number of neighbours yielding clearly more power on large scales, and less 
power due to noise on small scales - this is the expected effect of a better 
gradient accuracy due to a larger number of smoothing neighbours. How- 
ever, this advantage is quickly destroyed by the clumping instability. The 
bottom panel shows a histogram of SPH particle clump sizes determined 
with a friends-of-friends group finder, taking a linking length of 0.05 of the 
mean particle spacing. We see that the larger number of neighbours induces 
substantial clumping, reducing the number of independent sampling points 
and introducing large inhomogeneities in the kernel sums. 



In the top panel of Figure [8] we compare the velocity power 
spectra of these two simulations with the SI simulation. Note that 
at the resolution of 128 3 employed for these tests, the S2 run with 
512 neighbours is expected to have effectively the same mass- and 
spatial-resolution as the SI simulation with our default choice of 
64 smoothing neighbours. Interestingly, the power spectra look 
very similar on large scales, i.e. there is no noticeable improve- 
ment due to the higher number of smoothing neighbours at a fixed 
mass/spatial resolution. Only the small-scale noise is reduced when 
the number of neighbours is increased. However, if we look at early 
times after the driving has started (grey lines in Fig. [8}, larger dif- 
ferences are seen and the runs with a larger number of neighbours 
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Figure 9. Dependence of SPH turbulence results on the viscosity parame- 
terization. The top panel illustrates the effect of systematically varying the 
SPH viscosity strength. For lower a, the velocity power on large scales 
goes up, but the shape of the power spectrum does not improve. Note how- 
ever that this also increases the small scale velocity noise. However, if the 
viscosity strength is chosen too low, the power on large scales goes down 
again and the power spectrum is entirely noise-dominated in this case. The 
bottom panel compares different artificial viscosity schemes. Enabling the 
Balsara viscosity suppression factor improves the power spectrum, yielding 
a result similar to our a = 0.1 run. A further improvement is achieved if 
a time dependent vis cosity p arameterization is applied. For comparison, we 
include the results of Price] j20 12 1 for his 128 3 and 256 3 runs. The former 
is very close to our time-variable artificial viscosity run while the latter is 
comparable to our run with a = 0.01. The thin grey lines in both panels 
show the slope of the expected Kolmogorov power spectrum. 



do show more large-scale power, as expected for an improved ac- 
curacy. 

The bottom panel of Fig. [8]highlights the reason why this ad- 
vantage soon vanishes. Here we determine the size spectrum of par- 
ticle clumps formed in the runs at the end of the simulated time by 
applying a standard FOF algorithm with a small linking length of 
0.05 times the mean particle separation. Whereas in the S2 simu- 
lation with the default neighbour number essentially all particles 
stay isolated and no groups are found, this is very different in the 
runs with enlarged neighbour number. In the case of 512 smooth- 
ing neighbours, less than half of the particles remain isolated, with 
a large number of clumps containing multiple particles, up to ~ 10. 
This is the well-known clumping instability that frustrates attempts 
to beat down noise in the kernel sums by simply using a large num- 



ber of smoothing neighbours. Instead, the forming clumps reduce 
the effective resolution and the accuracy of the kernel interpolants. 

We now turn to the artificial viscosity parameterization, which 
is another area where one may hope that simple changes could lead 
to significant improvements in the results obtained for turbulence. 
In particular, the problematic damping of the injected turbulent en- 
ergy already on large scales hints that a reduction of the viscosity 
may help. A lower viscosity seems also warranted because in our 
subsonic regime shocks are not really expected, suggesting that ar- 
tificial viscosity may perhaps not be needed at all, or only at a mini- 
mal level. We have hence first repeated our default simulations with 
a serious of reduced settings of the artificial viscosity parameter, 
trying a = 0.1, a = 0.01, a = 0.001, and a = 0.0001. The top 
panel of Figure [9] shows the resulting velocity power spectra at S3 
resolution. Compared to our default S3 run, the large-scale power 
clearly increases when the viscosity strength is reduced, but at the 
same time the small-scale noise also drastically increased. In fact, 
we find that the energy dissipated in this noise-dominated regime 
is essentially invariant when the viscosity is varied (see Fig. [ 7 ]. 
While a larger artificial viscosity reduces the amplitude of the ve- 
locity noise, it also implies stronger viscous forces, such that the 
average work done against the viscous forces varies little. 

Eventually, however, the improvement of the large-scale re- 
sults when lowering the viscosity ends. Instead the results deteri- 
orate again when the viscosity is lowered to a = 0.001, or even 
a = 0.0001. In fact, for the latter case the large-scale result is even 
worse than for a = 1.0 whereas the small-scale noise is orders 
of magnitude larger. In this series of runs, it could be argued that 
a = 0.01 is “best” in the sense that the expected Kolmogorov slope 
of the velocity power spectrum is approximately reproduced over 
the widest range of scales. The large noise on small scales and the 
anemic ‘dip’ in the power spectrum at k ~ 200, which falls short 
of the expected shape of the dissipation range (see below), raise 
significant concerns though that the statistical properties of these 
turbulent motions suffer from significant noise contaminations. 

In the bottom panel of Fig. [9] we have instead enabled the so- 
called Balsara reduction factor for the viscosity in the presence of 
strong shear, and we consider a time-dependent viscosity parame- 
terization. Again, both of these changes lead to an increase of the 
power on large scales and more noise on small scales compared to 
S3, as expected due to the reduced viscosity. It turns out that the 
Balsara switch (S3-balsara) happens to lead to an extremely sim- 
ilar reduction of the effective viscosity as given by our run with 
a = 0.1 shown in the upper panel, but this is just by accident. 
Compared to these two runs, the simulations with time-variable vis- 
cosity give slightly more power on all scales, and are hence still less 
viscous overall. Incidentally, the run with out a shear viscosity lim- 
iter (S3-tav) agrees well with the result of Price 1 2012) for his 128 3 
run with time- variable viscosity (shown also in the figure, for com- 
parison), except on small scales due to Price’s choice of measuring 
the power spectrum of th e SPH-smoothed velocity field. However, 
the 256 3 result of Price ( 2012 1 with matching resolution lies still 



a bit higher, similar to our S3 result with a = 0.01. This is pre- 
sumably caused my small differences in how the time-dependent 
viscosity is parameterized. 



3.6 Reynolds numbers 

In a response paper to our original submission the present work, 
|Price| ( [2012j ) argued that our SPH results for subsonic turbulence 
can be naturally explained simply by our artificial viscosity set- 
tings. In fact, he suggests that SPH effectively yields a solution of 
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Figure 10. Velocity power spectra in SPH for different numerical settings 
but equal Reynolds numbers. In the top panel, we compare results for dif- 
ferent numerical resolutions, ranging from 64 3 (SI) to 256 3 (S3), and dif- 
ferent artificial viscosity settings, such that the Reynolds number estimated 
according to equation CD is the same and a very similar result would in 
principle be expected. In the bottom panel, we show instead results at a fixed 
resolution of 128 3 (S2), but here the Mach number is varied and the viscos- 
ity parameter is adjusted such that the estimated Reynolds number stays 
the same. An approximate self-similarity of the spectral shape is at best 
obtained only for large Mach numbers. For comparison, we also show cor- 
responding results for different Mach numbers for the moving-mesh code. 
Here the dynamic range of the inertial range is to a good approximation set 
by the grid resolution and is invariant with the Mach number. 



the Navier-Stokes equation with a kinematic viscosity v that can be 
estimated as 



V ~ a Vsia-h , 

10 8 



(16) 



where a is the artificial viscosity coefficient, v 3 i g ~ c a is the signal 
velocity, and h is the SPH smoothing length. He further argues that 
the Reynolds number of the SPH turbulence simulations is then 
given by 




Figure 11. Shape of the dissipation range in a moving-mesh simulation 
with physical viscosity, and in the SPH subsonic turbulence simulations of 
[Pncel (2012) . The dashed lines show fits to the velocity power spectrum 
based on equation m corresponding to Reynolds numbers 2100, 1000 
and 540 for AREPO and the 256 3 and 128 3 SPH results, respectively. The 
green dotted line is the expected shape for TZ e = 2100 turbulence, shifted 
in amplitude to fit part of the SPH result with 256 3 particles. It is clear that 
the dissipation range of the SPH results is not consistent in detail with a 
Kolmogorov cascade. We note that |Price](2012) quotes a Reynolds number 
of 6000 for his 256 3 result. 



is reduced away from shocks, and that our argument that gradient 
errors are responsible for a flawed Kolmogorov cascade would be 
“incorrect”. 

Actually, it is these assertions that are incorrect, as we show 
now. If the Reynolds number of SPH was indeed simply given by 
equation (17) as suggested by |Price| (2012) , we would expect in- 
variant results (in the sense of Reynolds-number similarity) for 
the turbulent power spectrum in different simulations when the 
Reynolds number and the driving are kept constant. This is how- 
ever not the case, as is shown in Figure[l0] Here we compare in the 
top panel a series of different SPH simulations in which the reso- 
lution is progressively increased by factors of 2, and the viscosity 
parameter a is correspondingly changed such that the (estimated) 
kinematic viscosity and Reynolds number stay fixed in all cases. 
Nevertheless the numerical results for SPH do not line up. In the 
bottom panel of Fig. [TO] we vary the Mach number of the turbu- 
lence at fixed resolution, again adopting different viscosity settings 
such that the Reynolds number according to equation (17) should 
stay the same. Again, an accurate universality of the shape of the 
SPH turbulence is not recovered, although there is a hint that this 
may work better for large Mach numbers. 

Another important point to make is that not only the inertial 
range of Kolmogorov turbulence for a Navier-Stokes flow is uni- 
versal, the dissipation range is universal as well (e.g. |Pope|2000) . 
In fact, experiments demonstrate that the energy spectrum can be 
well described by 



TZe = 



LV 

v 



10 L 
a h 



M. 



(17) 



Adopting as characteristic velocity scale the rmv-velocity of the 
turbulent velocity field (i.e. V = Mc s ), and for the length scale L 
the box size, |Price|(2012| ) estimates a Reynolds numbers of lZ e — 
6000 for his 256 3 run. He also claims that “turbulent flow with 
a Kolmogorov spectrum can be easily recovered” in SPH when a 



E(k) = Ce 2/3 k~ 5/3 f v (kri), 



(18) 



where e is the dissipation rate, and 



«/ 3 ' 1/4 
^ 1 - 



(19) 



is the Kolmogorov scale. The function f v (x) is universal and well 
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fit by 



f v (x) = exp — /3[(at 4 + c 4 ) 1/4 - ej) , 



( 20 ) 



with c ~ 0.4 and /3 ~ 5.2, and the value of the Kolmogorov con- 
stant, C ~ 0.5, is universal as well. Note that this means that 
for quasi-stationary turbulence, where the energy injection rate is 
equal to the dissipation rate, not only the inertial range is specified 
but also the shape of the spectrum in the dissipation range is fully 
specified if the kinematic viscosity is known. 

In Figure[TT] we ploMhe velocity power spectra for the 256 3 
and 128 3 SPH runs of 



Price 



2012 1 , and we compare them to a run 



with our Navier-Stokes solver in AREPO (Munoz et al. 2012), with 
v = 0.00015 at a resolution of 256 3 cells. The expected Reynolds 
number of this mesh-based simulation with prescribed physical vis- 
cosity (physical + numerical viscosity) is 7Z e ~ 2100. By plotting 
the spectrum using a linear scale for k, the dissipation range is em- 
phasized and becomes approximately a straight line. We readily see 
that AREPO provides an excellent fit to the expected shape (dashed 
line), as described by equations ( | 1 8| ) and ( |20| ), only a small devi- 
ation due to the bottleneck bump is present. However, neither the 



128 3 nor the 256 3 SPH results of Price 



(2012' provide a reasonable 



fit to the shape expected for Kolmogorov turbulence in the dissipa- 
tion range. This invalidates the claim that the subsonic turbulence 
results of |Price| ( |2012| ) (which are quite consistent with our own 
for the similar viscosity settings, see Fig. [9]i are consistent with 
Kolmogorov turbulence. It also shows that the Reynolds numbers 
quoted by |Price| < [20 1 2j ) are incorrect and too high by a factor of 
~ 6; in fact, his 256 3 run has at most 1Z E ~ 1000, and the 128 3 run 
about half that, but since the results do not accurately correspond to 
the expected Navier-Stokes solutions for these Reynolds numbers, 
these values have to be taken in any case with a grain of salt. We 
also want to remark that the dynamic range expected for the iner- 
tial range of Kolmogorov turbulence is of order p/L ~ 1Z E 3 ^ . A 
Reynolds number of 6000 should hence allow up to p/L ~ 680 - 
clearly infeasible with 256 points per dimension. 

Finally, we turn to the issue of gradient errors in SPH, which 
according to |Price| ( |2012| l play no role in the turbulence results. To 
examine this point, we carry out a simple experiment of a decay- 
ing large scale solenoidal velocity field. We set up a simple ran- 
dom realization of a number of solenoidal velocity modes on large 
scales (the largest 70 fc-modes in the box), with a ex fc -5 ^ 3 energy 
spectrum and normalized such that the resulting rms-velocity cor- 
responds exactly to a A4 = 0.3 Mach number. We now compare 
the time evolution of this field in three different simulations, with- 
out applying any driving. We expect that the large-scale shearing 
motions in the initial conditions will transfer some of their energy 
to smaller scales with time, and that the damping of these motions 
by numerical viscosity effects will eventually dissipate the kinetic 
energy of the system. One of the three codes we try is our moving- 
mesh code AREPO. The second is standard SPH with a low viscos- 
ity setting of a = 0.01. The third simulation uses the very same 
SPH code, but with the only difference being that the pressure in 
the equation of motion is replaced by P — » P — Po, where the 
constant Po = poc 3 is taken to be the background pressure. In 
principle, such a constant pressure offset should not change any- 
thing, as pressure gradients are unaffected and hence no effect on 
the dynamics is expected. However, SPH’s equation of motion has 
a ‘zeroth-order" error in the pressure gradient which is proportional 
to the pressure itself (e.g. Quinlan et al.||2006| |Read et al.||20T()| 



Gaburov & Nitadori 2011 |Amicarelli et al.||201 1[ ). For example. 
Read et al. (2010|> show that the actual acceleration acting on a 



SPH particle corresponds to 



dv» Pi (ViVO Pi 

-j— - -- — E 0 ,i h U(h), 

at hpi pi 



( 21 ) 



where E o,; is a dimensionless error vector, and Vi is a dimension- 
less error matrix. Only for E o,i = 0 and Vi equal to the identity 
matrix the correct equation of motion would be obtained. For an 
irregular particle distribution, one obtains however Eqj ~ 0(h ), 
meaning that the zeroth-order error is not easily reduced with better 
resolution. However, by subtracting Po, the average pressure of a 
particle is made close to zero, and with this trick the magnitude of 
the zeroth-order error should be greatly reduced. 

Indeed, if we look at the time evolution of the power spec- 
trum of our decaying velocity field (Figure[l2|l, we see a substantial 
difference between the two flavours of SPH. The standard version 
of SPH builds up small-scale velocity noise much more quickly, 
and the power on large scales decays more rapidly as well. In con- 
trast, the version of SPH with a pressure-bias in the equation of 
state manages to track the AREPO result on large scales for a much 
longer time. This hence proves that the small-scale noise created 
by gradient errors drains kinetic energy from large scales, effec- 
tively short-circuiting the Kolmogorov cascade. Notice that the two 
flavours of SPH examined here have identical viscosity settings, 
i.e. their only difference lies in the gradient errors in the SPH equa- 
tion of motion. Unfortunately, in general applications of SPH in 
astrophysics the subtraction of a constant background pressure is 
not readily possible, so the problem of gradient errors is present 
irrespective of the artificial viscosity settings. 

We note that the problematic influence of noise and gradi- 
ent errors in SPH has also become apparent in recent tests of the 
Gresho vortex problem ( |Springel||20 1 0 Read & Hayfield 201 1| >. 



The stable vortex flow in this set-up ( Gresho & Chan| 1990 1 tends 



to decay relatively quickly in standard SPH, but |Read & Hayfield] 
( |201 1| ) recently showed that the error and the convergence rate can 
be greatly improved if a higher-order gradient estimate, based on a 
much larger number of SPH smoothing neighbours and a different 
kernel shape, is used. 

It is also interesting to consider the computational cost re- 
quired to reach a certain effective Reynolds number in SPH and 
in a mesh scheme. If we ignore the issue of gradient errors in SPH 
for the moment and assume that the Reynolds number is indeed 
given by equation CD- then the computational cost to reach a cer- 
tain Reynolds number scales as fcpu oc 1Z E . This is because the 
number of particles scales as /i~ 3 , and the number of required 
timesteps as hT 1 . In contrast, because in a mesh code such as 
AREPO the Kolmogorov dissipation scale is essentially given by 
the cell size, p ~ h, we expect lZ e oc ft -4 / 3 , implying that the 
required CPU-time to reach a certain Reynolds number scales only 
as fcpu oc 1Z/. In the limit of large Reynolds numbers this is a 
significant competitive advantage. 



3.7 Vorticity power spectrum 

Another interesting probe for turbulence is provided by the vortic- 
ity w = V x v of the velocity field. This is because vorticity is in 
principle a conserved fluid quantity in an ideal gas, where the vor- 
ticity field is locked into the gas similar to a flux-frozen magnetic 
field in ideal magnetohydrodynamics. New vorticity is introduced 
on large scales through our solenoidal driving field, and it is erased 
on small scales via dissipation, but vorticity production through the 
baroclinic term should essentially be absent in our flow due to our 
quasi-isothermal conditions. Analyzing the statistical properties of 
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Figure 12. Time evolution of the velocity power spectrum of a decaying large-scale solenoidal velocity field. Initially, only the ~ 70 largest modes are 
populated with random phases and an expected E(k) oc k~ 5 ' 3 energy spectrum, normalized such that the rms velocity corresponds to M = 0.3. We 
compare three different simulation techniques: AREPO (blue), ordinary SPH (green) with a low viscosity of a = 0.01, and “pressure-biased” SPH (red), 
where the only difference is that a constant pressure Po = poc 2 is subtracted in the equation of state. The resolution is 128 3 in all cases. We see that ordinary 
SPH builds up small-scale noise considerably more quickly than the pressure-biased version of SPH. The latter actually tracks the AREPO result on large and 
intermediate scales much longer than ordinary SPH. As the viscosity in both SPH versions is the same, the faster dissipation of the large-scale motions in 
ordinary SPH is due to the larger noise on the smallest scales, as the energy dissipation rate E^ ISB (k) oc k 2 E[k) is dominated by these scales. 



the vorticity field can hence provide complementary information 
about turbulence and the numerical properties of the employed sim- 
ulation technique. 



In Figure [13] we show our measurements for the vorticity 
power spectrum for the same simulations that we used in Fig. [5] 
for an analysis of the velocity power spectrum. Both our moving- 
mesh and fixed-mesh simulations show approximately a power-law 
rise of the vorticity with scale, until a rapid drop sets in at around 
the numerical dissipation scale. For fully developed isotropic tur- 
bulence we expect a power-law spectrum of the vorticity given by 
E w ~ k 2 E v (e.g. Zhu et al. 2010). In particular, Kolmogorov’s 



theory then suggests a slope that is 2 units shallower than that of 
the velocity spectrum, and hence rises towards smaller scales as 
E w oc k 1 ^ 2 . This expectation is borne out well by our measure- 
ments, providing an important consistency check. Also, the propor- 
tionality between E w and k 2 E v is resolved quite well, as shown by 
the dashed lines. 

However, as already expected based on the enstrophy maps 
shown in Fig. [4] the SPH results for the vorticity power spectrum 
are very different. In essence, only some large eddies resulting from 
the driving are present, apart from some vorticity power on small 
scales, which is likely in large part noise, and in any case is smaller 
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Figure 13. Vorticity power spectrum of SPH and AREPO, compared at a 
resolution of 256 3 . The thin grey line gives the slope expected for Kol- 
mogorov’s theory of incompressible turbulence, in which for fully devel- 
oped turbulence, the power spectrum of the vorticity is proportional to k 2 
times the velocity power spectrum. The AREPO result follows this expecta- 
tion very well on large scales, over the same range where also Kolmogorov’s 
velocity power spectrum is reproduced. In contrast, SPH shows a rapid fall 
of vorticity towards small scales; only on scales of order the SPH smooth- 
ing length a substantial vorticity bump is seen, but this is presumably largely 
due to the velocity noise on these scales. 




Figure 14. Volume-weighted density PDFs for subsonic turbulence. Shown 
are SPH. AREPO, and fixed-mesh simulations. The PDFs are averaged over 
5 snapshots taken at times t = 12.8, 16.0, 19.2, 22.4 and 25.6. 



than what is found in the mesh code. This hence corroborates our 
previous results, confirming that the mesh-based simulations do 
resolve subsonic turbulence with the expected physical properties 
whereas this is more problematic in the case of SPH. 

3.8 Density probability distribution function 

The final quantity we examine in our simulations of subsonic turbu- 
lence is the density distribution function. In Figure[l4] we compare 
the volume-weighted density PDFs for our A3, F3, S3, and S3-tav 
simulations. These distribution functions have been averaged over 
5 simulation snapshots evenly spaced in time between t = 12.8 
and f = 25.6 since the start of the simulations. 



We find that the mesh codes agree well in their density PDF. 
However, the AREPO code in moving-mesh mode shows slightly 
more high-density regions than when a fixed-mesh is used. This 
is consistent with our expectation that the moving mesh yields a 
slightly lower numerical viscosity and a better adaptivity to high 
density regions. Due to the subsonic conditions in these simula- 
tions, the differences are expected to be quite small though. The 
density PDF of the SPH simulation shows somewhat larger differ- 
ences. The comparatively viscous a = 1.0 result is slightly thinner 
and hence more strongly peaked towards the mean value than for 
the mesh-based results. This is despite the fact that noise in the 
SPH density estimates can be expected to broaden the intrinsic dis- 
tribution, but apparently this effect is negligible compared to the 
physical width of the distribution one expects for Ad ~ 0.3 turbu- 
lence. Interestingly, our SPH run with time-variable artificial vis- 
cosity produces a density PDF very close to the mesh-based results 
on the low-density side, whereas it falls short a bit in the tail of the 
high-density side. 



4 TRANSSONIC AND SUPERSONIC TURBULENCE 

Superficially, the results obtained thus far seem to be in conflict 
with previous reports that SPH can adequately represent highly su- 
personic isothermal turbulence. However, it is important to appre- 
ciate that the physical nature of supersonic isothermal turbulence 
is really quite different from the subsonic regime studied thus far, 
both in terms of the role of ram pressure versus thermal pressure, 
and in terms of the relevant dissipation mechanisms. In supersonic 
turbulence, we expect a network of shocks, which in the limiting 
case of fully kinetic turbulence is described by Burgers turbulence 
and not by the Kolmogorov theory. 

We here briefly examine how well our new moving-mesh code 
AREPO describes turbulence in the transsonic and highly super- 
sonic regimes, and whether the SPH results improve in this regime. 
The names and key parameters of the primary simulations we have 
carried out for these tests are listed in Table[3] A first impression is 
obtained by the maps in Figure [15] where we show slices through 
the density and kinetic energy density fields in our Ai = 8.4 sim- 
ulations using a moving mesh in AREPO, a fixed mesh, or SPH. 
Right away we notice a much greater similarity of the maps than 
found in the subsonic regime, with SPH apparently being able to 
resolve the turbulence in a way that is at least qualitatively similar 
to the mesh-based results. 

We examine this more quantitatively in Figure |T6] where we 
show the velocity and density power spectra for increasing Mach 
numbers, ranging from the transsonic to the supersonic regime. We 
include results for all three types of numerical simulations that we 
have carried out, comparing them always with an identical forcing 
field as a function of time. 

The simulations shown in Fig. |T6] were performed for Mach 
numbers Ai ~ 1.2 (top row), Ai ~ 3.5 (middle row) and 
A4 ~ 8.4 (bottom row). We clearly see that as the Mach num- 
ber increases, the SPH method does progressively better for the 
velocity power spectrum and in fact appears to eventually converge 
to the result obtained with the two mesh-based techniques. While 
for Ai ~ 1.2, there is still a very significant deficit of power in 
SPH except for the largest scales, this effect becomes significantly 
weaker for Ai ~ 3.5 and almost vanishes for Ai ~ 8.4. When 
one compares the velocity power spectra of the moving-mesh and 
the fixed-mesh calculations, one generally finds very good agree- 
ment on large scales but a noticeable difference in the small-scale 
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Figure 15. Visualizations of the state of the gas in slices through our M ~ 8.4 supersonic runs at 256 3 resolution. The top row shows the logarithm of the 
velocity field and the bottom row the logarithm of the density field. From left to right, simulations with a moving mesh in AREPO (A3-ml0), a fixed mesh 
(F3-ml0), and SPH (S3-ml0) are shown. 



behaviour. The dissipation scale of the moving-mesh code lies at 
slightly smaller scale, which can be interpreted as a signature of a 
higher effective resolution for the same number of fluid cells. This 
advantage most likely arises from the reduced advection errors in 
the moving-mesh approach. 

If instead the density power spectra are compared (right col- 
umn in Fig. 1 1 6| ), we find a qualitatively similar behaviour. For 
higher Mach number, the SPH result approaches that of the mesh- 
based simulations. There are however somewhat larger residual dif- 
ferences between all three techniques at the highest Mach number 
compared with the situation found for the velocity power spectrum. 
In particular, the density power spectrum for the moving-mesh code 
has a shallower slope and extends to higher k than both for the 
fixed-mesh and the SPH codes. We interpret this as a result of the 
better adaptive resolution of the moving-mesh technique. Direct fits 
to the power-law region of the density power spectra at M ~ 8.4 
return slopes of —0.54 and —0.43 for these 256 3 moving-mesh and 
SPH runs, respectively, clearly indicating a significant difference. 
However, we note that the shape of the density power spectrum is 
relatively sensitive to resolution, as we show next, so these slopes 
are not numerically converged. 

In Figure[T7] we show a resolution study for our moving-mesh 
and fixed-mesh simulations for the cases of Ai ~ 3.5 and Ai ~ 
8.4 turbulence, considering results both for the velocity and density 
power spectra. It is seen that the velocity power spectra agree nicely 
between the moving-mesh and the fixed-mesh code at large scales, 
but that the effective resolution of the moving-mesh code is higher 
at a given number of resolution elements. When the resolution is 



improved, the power-law region corresponding to the inertial range 
is extended towards smaller scales, without a significant change in 
slope. A small bottleneck effect still seems to be present, but at a 
much smaller level than in the subsonic regime. 

In contrast, convergence in the density power spectrum is 
more challenging, as seen in the right column panels of Fig. HD 
Here low resolution can easily lead to an overestimate of the slope 
due to a fairly prominent bottleneck effect. Interestingly, in the 
A4 = 8.4 case we see that the slope of the fixed-mesh F4-ml0 sim- 
ulation is accurately reproduced already by the A3-ml0 simulation, 
at an almost one order of magnitude smaller number of resolution 
elements, and similarly for the F3-ml0 and A2-ml0 pair of simu- 
lations. This can be attributed to the adaptive nature and the lower 
advection errors of the moving-mesh approach compared with the 
fixed-grid Eulerian method. 

It is also interesting to examine the energy dissipation power 
spectra of the different simulation techniques in the highly super- 
sonic regime. This is done for our highest Mach number run in 
Figure |T8] All the simulations show a relatively broad distribu- 
tion of dissipation as a function of scale, which is quite different 
in character compared to the narrower distribution encountered in 
the subsonic case. This can be interpreted as a result of the different 
physical nature of the dissipation in this supersonic regime, which 
occurs primarily through a complex network of shock surfaces, and 
is hence not restricted to a small range of length scales. It is also in- 
teresting to note that SPH and the moving-mesh code show quanti- 
tatively a quite similar result, whereas the fixed-mesh method gives 
higher dissipation on nearly all scales. We argue that this is due 
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Figure 16. Turbulence power spectra of SPH, AREPO, and a fixed-mesh code for approximately sonic (top panel), mildly supersonic (middle panel) and for 
highly supersonic driven isothermal turbulence (bottom panel). The left panels show the velocity power spectra, the right panels the density power spectra. 
The dashed lines show fitted power-laws for AREPO (blue) and SPH (green). 



to the significant bulk motion present in the system, inducing en- 
hanced dissipation through advection errors in the fixed-mesh code. 
This is simply not present in this form in the two Lagrangian meth- 
ods, which are both Galilean-invariant schemes. 

Finally, Figure[l9]takes a look at the volume-weighted density 
probability distribution function (PDF) in the high Mach number 
case. We compare the PDFs of moving-mesh, fixed-mesh and SPF1 
simulations at the 256 3 resolution. The shape of all three results is 
described reasonably well by a log-normal distribution. Flowever, 



the fixed-mesh simulation shows a higher probability at the low 
density end and has the largest width of the distribution for this rea- 
son. The SPF1 simulation tends to give higher probability at the high 
density end, which is a very similar behaviour as found in |Price &| 
Federrath (2010). The moving-mesh run has an overall very simi- 
lar distribution as the SPF1 run, except for being slightly wider. To 
the extent that a better representation of the high-density tail is ad- 
vantageous in science applications of supersonic turbulence (which 
can be argued is particularly true in studies of star formation), the 
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Figure 17. Resolution study of the supersonic velocity and density power spectra for our fixed-mesh and moving-mesh simulations. The top panels show 
the runs at Mach number A4 ~ 3.5, and the bottom ones at A4 ~ 8.4. At large and intermediate scales, the moving-mesh runs correspond roughly to the 
fixed-mesh runs at 2 3 times higher resolution. 





logp/po 



Figure 18. Dissipation power spectra for the supersonic runs A3-ml0, F3- 
mlO and S3-ml0 at A4 ~ 8.4. As in Fig.[7]the velocity power spectra are 
plotted as dashed lines. 



Figure 19. The volume-weighted logarithmic density PDFs for our highly 
supersonic runs at A4 ~ 8.4, as labeled. The PDF is averaged over two 
snapshots at times t = 0.2 and t = 0.3. 



moving-mesh technique hence works at least as well as SPH, and 
arguably better than a fixed-mesh technique. 



5 DISCUSSION AND CONCLUSIONS 

Perhaps the most important question prompted by our results is 
why SPH behaves so badly in the subsonic regime. |Price| |20i2| 
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argued that the culprit lies simply in the artificial viscosity param- 
eterization, and that schemes that dynamically reduce the viscos- 
ity away from shocks (e.g. |Morris|1997[|bolag et al.|2005[|Cullen| 
|& Deh nen 2010) do much better and have no problem to repro- 
duce a Kolmogorov cascade also in the subsonic regime. While we 
agree that excessive artificial viscosity can compromise the results 
of SPH, particularly in the subsonic regime where this will show 
up more readily, this is by no means the complete story. Instead, 
we have demonstrated in this study that gradient errors inherent 
in standard formulations of SPH (‘classic SPH’) do play a major 
role as well. They seed small-scale velocity noise in shear flows 
on all scales, and this noise is dissipated away by the viscosity of 
the scheme. Since SPH is an energy-conserving scheme, this ef- 
fectively short-circuits part of the energy transfer cascade in Kol- 
mogorov’s theory of turbulence. 

The concern that the large subsonic noise in SPH may cause 
substantial accuracy problems in the treatment of fluid instabilities 
has recently been emphasized by a number of authors (Springel 
|20T0||Read~et al.|2010||Abel|2011| >. Also, numerous studies have 
pointed out that the standard approach followed in SPH for con- 
structing derivatives of smoothed fluid quantities involves rather 
large error terms, especially for the comparatively irregular particle 
distributions in multi-dimensional simulations. One problem lies 
in a lack of consistency of the ordinary density estimates (which 
do not conserve volume, i.e. the sum of rrn/ pi is not guaranteed 
to add up to the total volume), and another in a low order of the 



gradient estimate itself (e.g. Quinlan et al. 2006, 


Read et al.|2010, 


Gaburov & Nitadori 2011; Amicarelli et al. 201 1 


). In practice, this 



means that there can be pressure forces on particles even though 
all individual pressure values of the particles are equal, a point em- 



phasized in a recent study by |Abel| ( |201 1[ >. But if this is the case, 
spurious jittering motions of particles can be readily triggered even 
for a vanishingly small large-scale pressure gradient. 

In order to demonstrate this point explicitly and quantify the 
typical noise in the pressure gradient estimates of SPH and AREPO, 
we have carried out a simple experiment. To this end we used the 
particle coordinates x t of the last snapshot of our S3 subsonic sim- 
ulation run, which is representative for the particle distribution typ- 
ically encountered in SPH in this situation. We then assigned en- 
tropies to the particles (taking their density estimate into account) 
such that the pressures Pi = P(xi) of individual particles were 
given by the analytic pressure profile 

P(x) = P 0 q x + r, (22) 

which is a simple linear gradient in the q-direction (our results are 
independent of the actual orientation of this vector) with a constant 
pressure offset r. The SPH estimate for the pressure gradient was 
then inferred from the particle acceleration asm computed by the 
SPH code as 



VP = — osph p, 



(23) 



which is the relevant quantity that ultimately enters the discretized 
equation of motion. We can then consider the relative error of these 
SPH pressure gradient estimates with respect to the known analytic 
gradient. We define the corresponding errors as 



Crel 



\VP-P 0 q\ 

\Poq\ 






■VP 

IqIIvpi “ 



(24) 



and show them as scatter plots for a random subset of the points in 
Figure[20| 

For comparison, we also carried out the equivalent procedure 
for the AREPO code, based on the same particle coordinates. The 
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Figure 20. The top panel shows a scatter plot of the relative errors of pres- 
sure gradient estimates in a simple test set-up. The bottom panel shows the 
corresponding errors in the direction of the estimated pressure gradients. 
The errors of AREPO are at the level of machine precision. However, SPH 
shows severe errors in the pressure gradient estimates, with a size that de- 
pends on the magnitude of the pressure itself. If a constant offset of r = 1 is 
added to the pressure profile, the gradient errors are about ten times larger, 
exceeding even 100%. 



resulting errors are also shown in Figure [20] AREPO clearly cal- 
culates the pressure gradients highly accurately, both in magnitude 
and angle. In fact, AREPO’s gradient estimate is second-order ac- 
curate, independent of the distribution of points (fSpringel| [2010| ), 
implying that a linear gradient should be reproduced essentially to 
machine precision, which we find is also the case here. In contrast, 
SPH shows a huge scatter in both error measures. In fact, the mag- 
nitude of the absolute error can in extreme cases be up to twice as 
large as the value of the gradient itself, and also the angular errors 
are significant. Furthermore, the errors rise with the magnitude of 
p. If a pressure offset r = 1 is added to the pressure profile, the rel- 
ative error is increased by about an order of magnitude, which can 
be understood in terms of a higher Eq error in equation ( |21^ . We 
have also repeated the calculation for increased numbers of neigh- 
bours, which reduces the error, albeit only weakly. We note that 
these large errors occur for a rather simple problem - a spatially 
constant gradient. This makes it clear that standard SPH has com- 
paratively low-order accuracy for smooth flow. 

We thus think that the problems of SPH in resolving subsonic 
turbulence are serious. It is unlikely that they can be solved by 
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just increasing the resolution, reducing the artificial viscosity, or 
the number of smoothing neighbours. This is because changing the 
artificial viscosity parameterization does not improve the gradient 
estimates, and will hence not be able to resolve the underlying prob- 
lem. Furthermore, obtaining better gradient estimates through a 
larger number of SPH smoothing neighbours is efficiently blocked 
by the clumping instability, at least for the ordinary kernel of clas- 
sic SPH. In any case, what appears to be needed for better results in 
the subsonic regime are better gradient estimates. Some extensions 
and improvements of the standard SPH formulation that go into 



this direction have already been proposed (e.g. |Price|2008||Hefi &| 



Springe! 


2010 


Cullen & Dehnen 2010; Read et al. 2010 Read & 


Hayfield 


2011 


Abel 2011). It will remain to be seen whether any 



of them provides a robust and generally applicable alternative to 
standard SPH. 



We should clarify that despite the large errors in gradient esti- 
mates, it remains true that SPH has very good conservative proper- 
ties. This feature allows it to still produce physically sensible fluid 
behaviour in many situations despite the subsonic noise, especially 
in shock-dominated regimes where the accuracy of gradient esti- 
mates is much less important. Our results hence justify the appli- 
cation of standard SPH in studies of supersonic turbulence, pro- 
vided the Mach number is really high. At the same time, our results 
raise significant concerns for applications of SPH in regimes where 
subsonic phenomena such as turbulence are important. This is for 
example expected to be the case in cosmological structure forma- 
tion. Indeed, recent studies have already presented evidence that the 
accuracy problems of SPH in the treatment of the generation of tur- 
bulence and of fluid instabilities such as the Kelvin-Helmholtz in- 
stability affect galaxy formation directly jVogelsberger et al.|201 1| 
|Keres et al.|201 l[|Sijacki et al.|201 f] >. 

Another important conclusion from our results is that the new 
moving-mesh code AREPO is highly competitive with the accuracy 
of ordinary Eulerian mesh-codes for studies of turbulence. In the 
subsonic regime it produces essentially equivalent results as ordi- 
nary Eulerian codes, with a slightly reduced dissipativeness of the 
scheme. In the supersonic regime it however features a higher ef- 
fective resolution at the same number of resolution elements. In 
particular, the velocity and density power spectra can be traced to 
smaller scales, and there is generally less dissipation as a function 
of scale due to reduced advection errors. If just pure hydrodynam- 
ics without self gravity is considered, a moving-mesh calculation 
with AREPO is however more costly than a corresponding fixed- 
mesh or SPH calculation with the same number of resolution ele- 
ments. It is clear that an accurate description of turbulent gas mo- 
tions is highly desirable for a versatile astrophysical code, and we 
have shown here that AREPO is able to meet these requirements. 
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